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We present a technique for calculating free-energy profiles for the nucleation of multicomponent structures 
that contain as many species as building blocks. We find that a key factor is the topology of the graph 
describing the connectivity of the target assembly. By considering the designed interactions separately from 
weaker, incidental interactions, our approach yields predictions for the equilibrium yield and nucleation bar¬ 
riers. These predictions are in good agreement with corresponding Monte Carlo simulations. We show that a 
few fundamental properties of the connectivity graph determine the most prominent features of the assembly 
thermodynamics. Surprisingly, we find that polydispersity in the strengths of the designed interactions stabi¬ 
lizes intermediate structures and can be used to sculpt the free-energy landscape for self-assembly. Finally, we 
demonstrate that weak incidental interactions can preclude assembly at equilibrium due to the combinatorial 
possibilities for incorrect association. 


Building nanostructures out of multiple, distinct com¬ 
ponents offers enormous possibilities for high-fidelity 
manufacturing at the molecular level. Such ‘addressable’ 
structures are fundamentally different from conventional 
crystals or ordered clusters, since every building block is 
distinct and thus occupies a specific location in the tar¬ 
get structure. Because the interactions between building 
blocks are specified independently, it is possible to design 
finite-sized, three-dimensional structures that assemble 
nearly error-freej^^— An experimental proof-of-principle 
can be found in self-assembling DNA ‘tiles,’ which use 
the hybridization of complementary DNA sequences to 
construct complex structures consisting of hundreds of 
subunits from a single soup of monomers.— Simulation 
results have shown that such one-pot self-assembly can 
succeed with highly simplified model subunits that lack 
the molecular details of DNA tiles, suggesting that simi¬ 
lar design strategies should be widely applicable.— 

In this communication, we present a simple and effi¬ 
cient method for predicting the free-energy landscape for 
the self-assembly of addressable structures directly from 
the graph describing the target structure. Our approach 
allows us to predict the assembly yield and nucleation 
barriers quantitatively for any multicomponent structure 
in which directional, designed interactions stabilize the 
target assembly. We demonstrate the accuracy of this 
method by comparing our predictions with lattice Monte 
Carlo (MC) simulations of a DNA-tile nanostructure.— 
The good agreement suggests that our approach can be 
directly applied to any experimental system with des- 
ignable interactions. 

Our method relies on the observation that the designed 
interactions in the target structure are typically much 
stronger than any incidental associations between sub¬ 
units that should not be connected in the final assem¬ 
bly. These designed interactions can be represented by a 
graph, G. The vertices of this graph are all unique, since 
each vertex corresponds to a particular subunit. The 
edges indicate designed interactions that stabilize the 
target assembly. An example 26-subunit structure and 
its associated connectivity graph are shown in Figure [1] 



FIG. 1. (a) An example three-dimensional DNA-tile struc¬ 
ture in which all 26 subunits are distinct. For comparison 
with simulations,— the monomers are constrained to a cubic 
lattice with the pegs oriented toward one of the four nearest- 
neighbor sites. Every designed interaction between adjacent 
subunits is distinct, although the subunits can associate in 
any of the three possible dihedral configurations, (b) The 
connectivity graph, G, representing the designed interactions 
between arbitrarily labeled subunits, (c) A subgraph, g, of 
the target structure and its topological properties, as defined 
in the text. Bridges are highlighted in orange, and edges adja¬ 
cent to g are shown in gray, (d) The logarithm of the ‘density 
of states’ of fragments with E edges and V vertices. The ver¬ 
tical axis indicates the number of linearly independent cycles 
in fragments within each set. 


Initially, we consider only designed interactions that are 
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correct for error-free assembly of the target structure. 
Later, we shall discuss how to account for incidental in¬ 
teractions. 

In order to describe the assembly of the target struc¬ 
ture, we must consider the relative stability of every pos¬ 
sible correctly bonded partial structure (‘fragment’). In 
terms of our graph, a fragment corresponds to a con¬ 
nected subgraph of G. Fragments that differ only in the 
labels of their vertices are distinct because the vertices of 
G are all unique. Even though the total number of frag¬ 
ments grows exponentially with the number of edges in 
the target structure, limiting our attention to correctly 
bonded clusters ensures that the set of fragments is finite. 

We can deal with the enormous number of partial 
structures by grouping the fragments, {g}, into sets with 
the same number of edges, E{g), and vertices, V{g). 
While all fragments in a set h{E, V) may not share the 
same topology, they are likely to have similar thermody¬ 
namic properties. We then count the number of frag¬ 
ments in each set, \h{E,V)\, statistically by applying 
the Wang-Landau flat-histogram algorithm^, to the state 
space of connected subgraphs of G. In this algorithm, 
a stochastic trajectory hops among graphs in the frag¬ 
ment state space, making transitions between fragments 
that differ by a single edge. For every visited fragment g, 
we calculate the number of ‘bridges’ (edges that, if cut, 
would break g into two disconnected subgraphs), B{g), 
and ‘leaves’ (bridges that attach a single vertex to the 
rest of g), L{g) (Figured]:). The removal of a bridge that 
is not also a leaf would disconnect g into two nontrivial 
graphs and is thus not allowed. We also find all edges 
that are adjacent to 5 in G and denote the number of 
such edges by A{g). In order to obey detailed balance 
in this state space, randomly proposed edge additions or 
removals to transition between graphs g and g' are ac¬ 
cepted with probability 


Padd {g g') = min 1 , 

remove 


\HE{g),V{g))\ n±(g) ~ 
\h{E{g'),V{g'))\n^{g')\ ’ 


where n+{g) = A{g) and n_(g) = E{g) - B{g) -f L{g). 
When leaves are added or removed, the accompanying 
free vertex is attached or discarded from g' as well. By 
updating \h{E, V)\ according to the algorithm described 
in Ref. [tI, we can calculate the number of fragments in 
each set efficiently and to arbitrary precision. 

The resulting ‘density of states’ is an intrinsic prop¬ 
erty of the target structure that determines the most 
important features of its assembly at equilibrium. For 
reasons that will soon become clear, we have organized 
the {E, V}-sets in Figure [TJl according to the number 
of edges and the number of linearly independent cycles, 
G = E — V + 1. Conveniently, \h{E, V)| need only be 
calculated once for a particular target structure. 

In a dilute solution with many copies of each com¬ 
ponent, we can treat the mixture as an ideal solution 
of fragments. It is expedient to work in the grand- 
canonical ensemble, in which monomers of each compo¬ 
nent can be exchanged with an infinite reservoir at con¬ 


stant chemical potential. For simplicity, we assume that 
these chemical potentials are chosen such that each type 
of monomer is present in the same concentration, p, but 
in general, these concentrations can be different. Because 
all fragments {p} are in chemical equilibrium with the 
monomer reservoirs, the dimensionless grand potential is 
simply the sum of the fugacities {zg} of all fragments 
{ 5 }: — InS = —J2g^g (see, e.g. Ref. @). The sum over 
fragments can then be replaced by a sum over sets of 
fragments /i(i?, V), where each term is accompanied by 
the average fugacity of the fragments in the set, ze,v' 

Zid=lnS = ^|M£^,I/)|z^5,y. (2) 

E,V 

For every fragment, the associated fugacity can be de¬ 
termined directly from the fragment graph. First, there 
is an attractive contribution due to the designed bind¬ 
ing energies: -^Y.bee(g) where (3 = l/ksT, £{g) is 
the edge set of g and {£&} are the bond energies. Sec¬ 
ond, there is an entropic cost to bring monomers into 
contact, Vlnp. Finally, there is an entropic penalty due 
to the loss of rotational entropy of bonded subunits. Ig¬ 
noring excluded volume interactions, the formation of a 
bond forces the interaction sites on the monomers to face 
one another, reducing the dimensionless entropy by Inge, 
where gc is the coordination number of the lattice. Each 
bond that does not correspond to a bridge costs a fur¬ 
ther In gd of entropy, where gd is the number of dihedral 
angles possible for two bonded monomers on the lattice. 
The average fugacity of fragments in the set h{E, V) is 
thus 


lnzE,v = El3eE,v + Vlnp (3) 

-{V - l)lngc -{V - Be,v - l)lngd, 

where /3e is the absolute value of the mean bind¬ 
ing energy. The dimensionless dihedral entropy 

Be y In gd = In ( g?^®^) is also an intrinsic prop- 

’ \ / geh{E,V) 

erty of G and is easily computed with a stochastic calcu¬ 
lation in the fragment state space.— 

We first consider the case where all designed bond en¬ 
ergies have equal magnitude, i.e. /3e_E,y = f3e'iE,V. We 
define the yield, rj, to be the grand-canonical average of 
the fraction of all clusters in solution that match the tar¬ 
get structure, 


_ (No) _ ZG 

~ Eg{Ng) ~ 


( 4 ) 


where Ng is the number of copies of fragment g. Our 
prediction for the equilibrium yield of the structure de¬ 
scribed in Figure [T] is shown in Figure [2],. Quite strik¬ 
ingly, the transition from zero to nearly 100 % yield over 
approximately 2 k^T suggests highly cooperative assem¬ 
bly. Adjusting the monomer concentration simply shifts 
the yield curve, with the 50%-yield bond strength, /Seso, 
changing linearly with Inp. Increasing the bond strength 
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FIG. 2. (a) Comparison of the predicted equilibrium yield of the example structure and the largest stable cluster size, Vmax, 
obtained from MC simulations with a single copy of each subunit. All monomers are present in the same concentration, 
p = 62“®. (b) Comparison of the predicted and simulated free-energy profiles with /3e = 12.07. (c) The statistical weights in 
Eq. (HI) at /3eso- Details of the MC simulation methods used in (a) and (b) are provided in Ref. [ 3 . 


beyond /3e ~ 14 results in perfect assembly because only 
designed interactions are considered at this point. The 
predicted yield curve coincides remarkably well with the 
largest stable cluster observed in MC simulations. 

The correspondence between the intrinsic properties of 
the connectivity graph and the equilibrium self-assembly 
of the target structure is immediately apparent from 
the free-energy profile shown in Figure Hb- Because we 
are interested in the progress toward complete assem¬ 
bly starting from any subset of components, the relevant 
free energy is a sum over all fragments consisting of V 
monomers, 

/3F(V) = -lnJ2lh(E,V)lzE,v. (5) 

E 

Exactly seven free-energy barriers, each corresponding 
to a linearly independent cycle in G, must be crossed 
in order to assemble the example structure from free 
monomers in solution. Again, we find good agreement 
with the results of MC simulations^^ The ascent on ap¬ 
proach to each barrier results from the recruitment of 
additional subunits by single bonds. Each steep descent 
corresponds to the completion of a new cycle, in which 
the additional bond compensates for the loss of rotational 
and translation entropy. The critical nucleus is V = 5, 
since six subunits are required to form a cycle in the ex¬ 
ample structure. Beyond the critical nucleus, the num¬ 
ber of monomers required to complete additional cycles 
can be determined directly from \h{E,V)\. Unsurpris¬ 
ingly, fragments with the greatest number of bonds per 
monomer are most stable and thus most likely to form at 
equilibrium, as shown in Figure HJ:. 

But what if the strengths of all designed interactions 
are not identical? If the designed interactions are instead 
chosen randomly from a distribution with mean —jSe and 
a finite variance a^, then we should expect the thermo¬ 
dynamic properties of the self-assembling system to fluc¬ 
tuate as well. Because the bond energies do not change 
during assembly, we must treat the {eb}’s as quenched 
random variables when taking thermal averages. The 
mean bond energy of the target structure is clearly self¬ 


averaging, since e{E{G),V{G)) —>■ e by the central limit 
theorem. In the case of partial structures, however, the 
effective mean bond energy fdcEy depends on the number 
of fragments in the set h{E,V). Because self-averaging 
is only meaningful in the context of free energies that are 
extensive in E{G)j^ the effective mean bond energy is 

/ScB.y = ( ^ In / exp ^/3eh \ \ . (6) 

\ \ b€Sig) 

The inner (■) is a thermal average over all frag¬ 
ments with fixed bond energies corresponding to 
each edge in G. Assuming a Gaussian distribu¬ 
tion of independently chosen {/?e&}’s, we would find 
ln(-)g/£; ^ e-f ^a^{l - [E - 1]/[E(G) - I]} if all com¬ 
binations of edges were allowed in the fragments. The 
requirement for subgraphs to be connected, however, in¬ 
troduces additional correlations among the bond energies 
in h{E,V). The outer (•) averages over the quenched 
random interactions, which introduces an additional cor¬ 
rection due to the finite number of bonds. The effec¬ 
tive mean bond energies of the example structure, with 
= 1 and tr^ = 2.5, are shown in Figure [3^. 

The physical consequence of this analysis is that poly- 
dispersity in the designed interaction energies tends to 
stabilize partially formed structures. This effect origi¬ 
nates from the exponential weighting of the bond ener¬ 
gies, which is most significant when the fragments are 
small and the number of fragments in a set is large. Be¬ 
cause relatively few fragments with a given number of 
edges have many cycles, energy polydispersity also tends 
to stabilize ‘floppy’ fragments. As a result, the yield 
curves broaden and shift to higher fie as the variance 
in the interaction energy distribution is increased (Fig¬ 
ure [3]:). Energy polydispersity also rounds off the peaks 
of the free-energy barriers, diminishing the signature of 
the cycle-dependence for very broad interaction energy 
distributions (Figures |3Ji-e). 

Finally, we can estimate the thermodynamic conse¬ 
quences of incidental interactions between subunits. In 
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FIG. 3. (a) Effective bond energies for each set of fragments, j3iE,v, averaged over Gaussian-distributed designed interactions 
with variance cr^. (b) The ratio of the contributions to the grand potential due to incidental interactions, 2in, and designed 
interactions, Zid. (b,inset) The statistical weights of incorrectly formed dimers of fragments with V and V' vertices at /3e5o 
for the case of identical bond energies, (c) Predicted average yields given polydisperse energies, with and without incidental 
interactions, (d) Predicted average free-energy profiles with polydisperse energies, (e) Gomparison of MC simulations with 
identical bond energies and a particular set of quenched energies. 


terms of the connectivity graph, incidental interactions 
occur between vertices in two fragments that both have 
adjacent edges but do not share an edge in G. Assuming 
that the mean incidental interaction strength, (iw, is sig¬ 
nificantly weaker than /3e, we can use a high-temperature 
expansion to estimate the effect of incidental interactions 
in a reference state of ideal fragments^ The additional 
possibilities for binding result in a decreased dimension¬ 
less grand potential, — InS = — 0(w^). The 

second term in this expansion accounts for all possible in¬ 
cidental interactions between fragments in the ideal ref¬ 
erence stater^ 


min(Ag,A ,) 


Zin — ^ ZgZg 

g,a' 


E 

A 


A, 


Ag> 


■) — 15 „ 


A!- 




'' ^p/3-w 

.2A-1 ’ 


where the sum over A accounts for the many in¬ 
cidental binding opportunities between the ‘multiva¬ 
lent’ fragments*^ Zjn can be rewritten as a sum over 
{E,V;E',V'}, with Ae,v = {Ag)g^h{E,v) determined 
from a stochastic calculation in the fragment state space. 

For successful assembly to occur at equilibrium, inci¬ 
dental interactions must be less probable than designed 
interactions, such that < Zid- Figure HJo shows that 
Zin/Zjd increases sharply as soon as multimeric frag¬ 
ments become populated (,5e ~ 11.3). As a result, the 


yield quickly drops to zero, even with extraordinarily 
weak incidental interactions (Figure [31:). Although the 
point at which Z^ exceeds Zjd is relatively insensitive to 
both j3w and the equilibrium assembly window nar¬ 
rows rapidly since larger clusters, which are numerous 
under conditions of high yield, are also the most suscep¬ 
tible to aggregation (Figure |3l3,mset). Incidental inter¬ 
actions therefore present a fundamental thermodynamic 
constraint for successful self-assembly at equilibrium. 

The method presented here is generally applicable to 
any addressable structure that is stabilized by specific, 
directional interactions. By considering the complete set 
of on-pathway assembly intermediates, this approach re¬ 
veals how the fundamental topological properties of the 
connectivity graph determine the free-energy landscape 
for self-assembly. This method is also spectacularly ef¬ 
ficient: the calculations required to generate Figure [^Id, 
for example, were six orders of magnitude faster than the 
corresponding MC simulations, and the fragment den¬ 
sity of states can be reused to compute free-energy pro¬ 
files at any temperature and concentration. Futhermore, 
generalizations of this method may be used to study 
nucleation in complex crystals and the assembly path¬ 
ways of heteropolymers with designed native structures. 
This graph-based approach therefore points the way to¬ 
ward improved designs for a broad class of self-assembling 
nanostructures. 
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